29th Monitoring Research Review: Ground-Based Nuclear Explosion Monitoring Technologies 


ADVANCED WAVEFORM SIMULATION FOR SEISMIC MONITORING EVENTS 

Don V. Helmberger 1 , Jeroen Tromp 1 , and Arthur J. Rodgers 2 

California Institute of Technology 1 and Lawrence Livermore National Laboratory 2 

Sponsored by National Nuclear Security Administration 
Office of Nonproliferation Research and Development 
Office of Defense Nuclear Nonproliferation 

Contract No. DE-FC52-06NA27319 1 ’ 2 


ABSTRACT 


Comprehensive nuclear-test-ban monitoring in terms of location and discrimination has progressed 
significantly in recent years. However, the characterization of sources and the estimation of low yields 
remains a particular challenge. As the recent Korean shot demonstrated, we can probably expect to have a 
small set of teleseismic, far-regional and high-frequency regional data to analyze in estimating the yield of 
an event. Since stacking helps to bring signals out of the noise, it becomes useful to conduct comparable 
analyses on neighboring events, earthquakes in this case. If these auxiliary events have accurate moments 
and source descriptions, we have a means of directly comparing effective source strengths. Although we 
will rely on modeling codes, ID, 2D, and 3D, we will also apply a broadband calibration procedure to use 
longer periods (P >5 s) of waveform data to calibrate short-period (P between 0.5 to 2 Hz) and 
high-frequency (P between 2 to 10 Hz) as path-specific station corrections from well-known regional 
sources. We have expanded our basic cut-and-paste (CAP) methodology to include not only timing shifts 
but also amplitude (f) corrections at recording sites. The name of this method was derived from source 
inversions that allow timing shifts between “waveform segments” (or cutting the seismogram up and 
re-assembling) to correct for crustal variation. For convenience, we will refer to these f-dependent 
refinements as CAP+ for (short period, SP) and CAP++ for still higher frequency. These methods allow the 
retrieval of source parameters using only P-waveforms where radiation patterns are obvious as 
demonstrated in this report and are well suited for explosion P-wave data. The method is easily extended to 
all distances because it uses Green’s function although there may be some changes required in t* to adjust 
for offsets between local vs. teleseismic distances. In short, we use a mixture of model-dependent and 
empirical corrections to tackle the path effects. Although we rely on the large TriNet array as a testbed for 
refining methods, we will present some preliminary results on Korea and Iran. 
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OBJECTIVE 


Our main objective is to develop methods for locating and identifying seismic events, focusing primarily on 
earthquakes but including mining blasts and explosions. To advance methodologies, we want to use our 
computational resources to address both deterministic approaches as well as to develop a better physical 
understanding of scattering where frequency-dependent coda magnitudes prove so effective. While the 
latter is apparently source-reduction pattern free (above some cut-off frequency), the direct phases are not 
and thus the combination should be even more effective in discrimination. In this report, we concentrate on 
the direct phases and produce a very well recorded set of events (2 < M < 5) which can be used for 
calibration and scattering models. 

RESEARCH ACCOMPLISHED 


The primary objective of this research is to use earthquake data to model and calibrate paths so that more 
accurate estimates of explosion strengths can be established for regions of interest. The basic idea of using 
events (M w > 4) to calibrate paths has been introduced earlier, i.e., Song and Helmberger (1998). They 
attempted to build on the empirical Green’s function method, Hartzell (1978), by generating pseudo 
Green’s functions that are less dependent on source excitation. The method uses the perturbation of 
individual generalized ray responses calculated from a ID model where each layer is divided into blocks 
and allowed to vary in velocity. The small shifts produced by these tomographic changes are assembled to 
force ID Green’s functions to fit data better using a simulated annealing algorithm. Application to a set of 
the Landers’ earthquake aftershocks is given in Helmberger et al. (2001) where the most variable velocities 
occur along the surface or just above the Moho. 

Several automated waveform modeling techniques are currently available for estimating source parameters 
based on fitting long-period (LP) data (periods > 5 s). Most methods rely on a library of ID Green’s 
functions which are compared against observed records to deduce moment, mechanism, and depth for 
M > 3.5. A method we prefer uses a cut-and-paste procedure (CAP) which allows P-waves to be fit 
independent of the stronger surface waves, where timing shifts are allowed between phases. This approach 
is somewhat simpler than those discussed above. These shifts can be used to derive a new 3D model as 
demonstrated in this report or stored as a tomographic-type map to be used in extracting parameters for 
smaller events. An example of the CAP analysis is given in Figures 1 and 2 for a moderate magnitude event 
occurring in Iran. Note that the Pnl portions of the vertical and radial segments have been enhanced on the 
left and fit separately. The numbers below each trace indicates the cross-correlations (top) and shifts 
(lower). The header line states the model, the strike (59°), the dip (80°), and the slip direction (40°). The 
depth is estimated to be 4 km, and has a magnitude M w of 5.2. The number after the station name indicates 
the distance in km. This event is well-recorded globally and becomes a natural starting place for using a 
combination of regional and teleseismic data to study events. 

A collection of this type of data is being assembled (Rodgers et al., 1999). However, we concentrate on the 
technique development using the Southern California Seismic Network as a testbed. We address five issues : 
(1) refining models and CAP validation, (2) two-station solutions, (3) extension of CAP to contain 
amplitude corrections, (4) directivity and cluster analysis, and (5) site characterization. 

1. Refining Velocity Models and Validation 

Using the CAP method, we have avoided propagational complexities by allowing shifts in record segments. 
Here we will use these shifts to refine crustal structure. The procedure for refining the upper crust is 
displayed in Figure 3. We have obtained the shift required to align the ID synthetics with the observed 
waveforms for about 200 events. The delays are given in terms of a spider diagram where the model 
produces early synthetics in general (red). The cross-correlations are displayed in similar fashion where 
most of the fits are good (blue). The delays for all of the events can be assembled and used to generate a 
tomography image (block-type) and slightly smoothed as displayed for Love waves. Earlier sensitivity 
studies indicate that the seismogenic layer containing the source controls the timing, Song et al. (1996). 
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Synthetics (upper right) can then be generated with a fast (2D) FD routine assuming the mechanisms 
obtained from the CAP solutions, Helmberger and Vidale (1988). The bottom seismogram was generated 
with FK-1D, Saikia (1994). Comparisons for the ID and 2D FD synthetics are also included for comparison 
with the data trace (top). These records are filtered to 5s which is the normal filter (LP) used in the CAP 
inversion. The spider diagram at the bottom (2D) indicates the improvement in timing of the new model 
although the fits remain unchanged. These calculations are fast, taking about 5 minutes per event on a 
desktop computer. Similar results are obtained for the other events and will be presented at the annual 
meeting. 

2. Locating and Modeling Regional Earthquakes with Two Stations 

We have introduced a new technique to retrieve source parameters including location by using established 
shifts, either from measurements or 3D models, Tan et al. (2006). The method conducts a grid search in a 
cube centered over the source target comparing synthetics with shifts to observed records to determine 
location, depth, origin time, and mechanism that best fits the 3-component seismograms. The method was 
tested on 28 events beneath the Tibetan Plateau against a Program for the Array Seismic Studies of the 
Continental Lithosphere (PASSCAL) array solutions. Here we test over 160 events against a large station 
distribution, with a sample of the average fits and one of the poorest (Figure 4). 

3. Extension of CAP to Shorter Periods 

The changeover from LP (> 5 s) to SP (< 0.5 s) waveform modeling is facing the inherent trade-offs 
between source complexity and structural heterogeneity. Under such circumstances, analyzing clustered 
events provides a practical way to “separate” the source from the structural effect. The feasibility is 
displayed in Figure 5, where we compare the records at the same station GSC from three clustered events of 
different sizes. The two smaller events are overwhelmed by noise in long-period (0.05-0.2 Hz) frequency 
band; however, they produce decent signals at higher frequencies, such as 0.5-2 Hz, or 2-8 Hz. In 
particular, at 0.5-2 Hz, the discrepancies between the records from the magnitude 4 (13938812) and 
magnitude 2 (13937632) events are mainly caused by their differences in size and focal mechanism, while 
the effects from the detailed rupture processes are considerably smaller. This suggests that short-period 
P-waves can be used for determining focal mechanisms of magnitude 2 events if we could calibrate the path 
effect using magnitude 4 events with known source mechanisms. The same is true when it comes to the 
magnitude 2 (13937632) and magnitude 1 (13939108) events at higher frequency (2-8 Hz). On the other 
hand, the smaller events can provide excellent empirical Green’s functions for studying the detailed rupture 
process of the magnitude 4 event. Such a “two way” calibration process is what we have demonstrated in 
our recent papers (Tan and Helmberger, 2007a and 2007b). 

The discrepancies between the observed P-waves and the synthetics are mainly manifested as amplitude 
differences. For a quantification purpose, we define the function of “amplitude amplification factor” (AAF) 
as 


AAF 


| u 2 ( t)dt 
J s 2 ( t)dt 


( 1 ) 


where u(t) and s(t) are the data and synthetics, respectively. The integration is over a 2 s window centered 
on the onset of the P-wave. It appears that the most anomalous AAFs occur for the stations in the basins. In 
particular, these stations are consistently characterized by large AAFs (> 1) on the vertical component, but 
small AAFs (< 1) on the radial component. This discrepancy between the vertical and radial components 
has been noted by many previous investigators (e.g., Savage and Helmberger, 2004) and is addressed later 
in this proposal. By comparing the AAFs derived from all the calibration events, particularly, the thrust and 
the strike-slip events, we found a large number (-70) of stations display stable and mechanism-independent 
AAFs, which suggests the amplitude discrepancy of the synthetic P-waves could be corrected to model the 
observations, hence determine the earthquake source parameters. 
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We invert the short-period P waves with a similar grid-search approach as in the long-period inversion, 
where we minimize the L2 norm of the misfit between the data and synthetics: 

e — | u(t) - AAF . s ( t ) ||. (2) 

The AAFs in equation (2) are taken as the averages of the AAFs derived from all the calibration events to 
account for the structural effect. The few stations that display rapid variations of the AAFs among the 
calibration events are discarded. The validation test with respect to the calibration events shows remarkable 
agreement between the results from short-period P-wave inversions and their LP solutions. The same 
strategy can be carried to higher-frequency bands, e.g., 2-8 Hz, and hence be applied to smaller events. To 
examine the AAFs at higher frequencies and their usefulness in determining smaller events’ focal 
mechanisms, we conducted an AAF calibration process at 2-8 Hz with 21 M w ~ 2.5 events from the Big 
Bear cluster. Although the site deviations are generally larger, the AAFs still appear stable and mechanism- 
independent. More importantly, when we conducted P wave inversions with these AAF corrections for the 
calibration events, their resulted mechanisms agree with those determined at 0.5-2 Hz solutions. A few 
examples of the resultant waveform fits are given in Figure 6. Note the nodal behavior of the few 
stations, such as NBS, HEC, BEL, BLA, RVR, and SVD. This preliminary study suggests that the 
P-wave radiation pattern is still preserved for a large number of stations at 2-8 Hz. As the AAFs for 
these stations are calibrated, they can be used for studying source mechanisms of magnitude 1 events. 
Thus, these small events can be used to study the directivity of events 2.5 < M w <3.5 using the same 
procedure given by Eq. (4). Note at these frequencies, we can study small explosions. 

4. Rupture Directivity and Cluster Analysis 

Most large events (M > 3.5) show considerable rupture directivity in terms of source duration variation 
with azimuth. Here we will take the magnitude 2 events as empirical Green’s functions (EGFs) to 
investigate the rupture processes of the 7 larger strike-slip events (M w > 3.5) from the Big Bear sequence. 
Rather than deconvolution (e.g., McGuire, 2004) or estimating finite-moment tensors (FMTs), we propose 
an alternative approach of forward modeling to retrieve the relative source-time function (RSTFs). Let d(t) 
and g(t) be the records from a large event (M w >3.5) and the associated EGF event at the same station, 
which can be related by the relative source time function, RSTF(t) of the large event as: 

d{t) = g(t)*RSTF(t). (3) 

Assuming a simple trapezoidal shape of RSTF(t) according to the ID Haskell model (Haskell, 1964), where 
a RSTF(t) can be parameterized as the convolution of two boxcars, featuring the rise time x\ and the rupture 
time x 2 , we can solve for the RSTF(t ) in a grid search manner by minimizing the misfit defined as 

e = jd (V) - AM 0 g(t) * RSTF (f)||, 
where RSTF (t) = T x (7) * r 2 (f). 

Here the parallels denote the L 2 norm. A M 0 is an amplitude scaling factor to account for the two events’ 
known difference in size and radiation patterns. We first use the event 13937492 to illustrate the whole 
process, which is also the smallest (M w ~ 3.5) among our analyzed events (Figure 7). 

5. Site Characterization 

Thus, it appears that the greatest barrier to modeling at short periods (< 2 s) or wavelengths less than 10 km 
is simply the scattering produced by the surface. This is no surprise to those researchers working with 
down-hole arrays, but it is surprising how well simple calibration works to correct for some of these 
complicated features. Although not discussed here, one of the greatest advantages of working with 
wavelengths less than a few km is that directivity effects are relatively obvious as displayed earlier. While 
the AAFs themselves are back-azimuth and distance-dependent, their ratios between the vertical 
and radial components appear constant and are possibly simple functions of site conditions 
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(Figure 8). Note the clustering of the AAF(v)/AAF(r) ratios at the various stations from events at 
different locations. There is also good correlation between the AAF(v)/AAF(r) ratios and the 
surface geology; the relationship between this ratio and the site condition will be pursued as more 
measurements are made. 

RECOMMENDATIONS 

In summary, we are adjusting our research efforts to address one of the primary issues facing the Air Force 
Technical Applications Center (AFTAC jmonitoring community, namely, discrimination at high frequency 
(HF) and yield estimation of small explosions. Since these events are P-wave rich, we have extended our 
regional modeling efforts to 8 Hz where the surface conditions at the recording station play a key role. To 
achieve this, we rely on empirical corrections of ID Green’s functions stepping downward in magnitudes as 
the fault-dimensions become shorter and shorter, avoiding rupture properties, i.e., 4’s are used to study 3’s, 
3’s to study 2’s, etc. This hybrid approach now appears to work for CAP (shifts only), CAP+ (corrected for 
amplitude, 2 to 0.5 Hz), and CAP++ (corrected for 2 to 8 Hz). These corrections are path-dependent and 
rely on earthquakes for calibration although the shifts can probably be obtained from tomography and site 
conditions from topography images. The next step is to apply these methods to explosions, including 
Nevada Test Site explosions, and compare with other methods, Coda magnitudes, etc. 
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Figure 1. GSN (red square), KZ 
and KN (red triangles) and IIEES 
(black triangles) stations in and 
around Iran. Blue dots indicate 
earthquakes from the National 
earthquake Information Center 
(NEIC) catalog (big dots: 
earthquake > M5, small dots: 

< M5). The red star is the 060507 
event modeled with the CAP 
algorithm for IIEES stations. 
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Figure 2. Application of CAP to the modeling of the 060507 event. This earthquake features a 

mostly strike-slip mechanism with some normal component. The focal depth is estimated 
to be around 4 km, as compared to 14 km from the IIEES catalog. The closest station is 
90 km away, so depth ca not be well resolved just from first arrivals. 
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Figure 3. The top panel displays the Love wave phase velocity perturbation in southern California 
with a selected source (star)-receiver (triangle) path, along which the 2D cross section is 
shown and the synthetics are compared against the data on the right. The lower panel 
summarizes the comparison between ID and 2D synthetics against the observed Love 
waves for a Big Bear event. 


86 







































































29th Monitoring Research Review: Ground-Based Nuclear Explosion Monitoring Technologies 


35°00' 


34°45' 


-iifi’srv 


-11 fi°1 fV 



Figure 4. Grid-search results of source parameters for two southern California events using the 

two stations GSC and PAS. The top panel displays the source-receiver geometry, and the 
lower panel displays the results. The white stars denote the best solution with the 1.2 
contour for the location uncertainty estimates. Also included for comparison are the 
locations from Lin and Shearer (2007) (red stars) and the mechanisms from the whole 
TriNet array (Tan, 2006). 
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Figure 5. Comparison of the GSC records of three clustered, but different-sized events from the 2003 
Big Bear sequences. The original broadband and filtered records of different frequency 
bands (0.02-0.2 Hz, 0.5-2 Hz, 2-8 Hz) are displayed from top to bottom. The short period 
comparisons are concentrated on the P-wave trains. 
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Figure 6. Comparison of the P-wave waveform fits from the 2-8 Hz short-period inversions for three 
events with different mechanisms. 
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Figure 7a 
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Figure 7b Figure 7c 

Figure 7. (a) The selected waveform fits (vertical P-waves) between the records from event 13937492 
(black) and the “synthetics”(red) constructed with EGFs from event 13937632. The RSTFs 
are given to the left. Plotted are the absolute amplitudes, except that a scaling factor of 1/4, 
1/2 and 2 has been applied to the stations JVA, PDU, and PLS, respectively, for the display 
purposes. The obtained best RSTFs for the stations are circled. Note the apparent 
azimuthal pattern of the RSTFs. 

(b) .The vertical P wave amplitude ratios between the event 13937492 (M ~ 3.5) and 
13937632 ( M ~ 2.5). The two events have similar strike slip focal mechanisms. The two 
gray lines display the strikes of their fault planes. Note the amplitude ratios for the 
southeastern stations are consistently larger than those for the northwestern stations, due 
to the rupture directivity effect of the bigger event. 

(c) The inferred rupture directivities from the Big Bear cluster. The arrows point to the 
rupture directions, while their lengths (for the solid ones) indicate the fault lengths. Note 
that the ruptured planes correlate well with the seismicity lineations, suggesting cross-over 
faults at depth. 
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Figure 8. (top) The comparison of the ratios between the AAFs for the vertical and radial 

components from six events vs. those from the 2003 Big Bear sequence. The stations are 
grouped and color-coded according to their site-conditions as shown in the bottom, 
(bottom) The source parameters of the events for comparison (from Tan (2006)) 
and the station locations are displayed on a geologically-determined 
site-condition map from Wills et al. (2000). Labeled is the representative V s 30 of 
the top 30 m for each category, and the expected range of V s . 
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